This notebook performs pair-wise comparisons of qPCR gene expression, normalized to GAPDH expression. It calculates delta Cq, delta delta Cq, and fold changes in expression. Additionally, it generates box plots (delta Cq), and bar plots (fold change expression).
Groups are named in the following fashion:
<life.stage>.<conditioning.treatment>.<acute.treatment>
This allows for parsing downstream.
NOTE: Below is the full set of groups for the entire experiment. For the current qPCR analysis, seed and spat do not have acute treatments; just conditioning treatments.
seed.control.ambient=c("29", "40", "55", "63", "69", "101", "119", "122", "155", "164", "187", "202", "209", "214", "233", "236", "275")
seed.control.high=c("42", "59", "60", "62", "86", "102", "140", "176", "177", "184", "192", "223", "234", "243", "244", "254", "264")
seed.treated.ambient=c("14", "48", "66", "72", "89", "115", "129", "138", "156", "182", "191", "201", "227", "239", "270", "277", "280")
seed.treated.high=c("15", "19", "24", "88", "92", "105", "111", "113", "120", "128", "161", "200", "211", "256", "257", "266", "285")
spat.control.ambient=c("11", "30", "36", "52", "77", "114", "134", "142", "144", "183", "193", "229", "230", "231", "240", "272", "287")
spat.control.high=c("27", "74", "93", "96", "97", "137", "143", "153", "168", "178", "189", "206", "262", "274", "282", "284", "289")
spat.treated.ambient=c("9", "13", "38", "46", "47", "121", "145", "151", "174", "194", "197", "198", "216", "235", "241", "252", "291")
spat.treated.high=c("6", "25", "50", "78", "124", "126", "131", "160", "163", "172", "220", "226", "242", "253", "296", "298")
juvenile.control.ambient=c("18", "57", "65", "75", "79", "104", "110", "123", "125", "171", "175", "205", "238", "273", "279", "293", "317")
juvenile.control.high=c("12", "39", "43", "49", "71", "130", "141", "146", "150", "170", "195", "297", "301", "324", "351", "355", "371")
juvenile.treated.ambient=c("1", "34", "64", "83", "98", "147", "152", "158", "162", "169", "188", "271", "295", "310", "357", "361", "381")
juvenile.treated.high=c("28", "53", "61", "73", "81", "106", "109", "139", "149", "173", "181", "213", "290", "302", "311", "364", "392")
adult.control.ambient=c("3", "5", "13*", "16", "17", "80", "87", "94", "148", "159", "179", "180", "250", "258", "268", "312", "326", "330", "334", "346", "360", "377", "379", "386")
adult.control.high=c("20", "23", "26", "32", "33", "67", "70", "90", "107", "132", "135", "157", "166", "186", "207", "215", "248", "316", "341", "344", "349", "382", "394", "395")
adult.treated.ambient=c("7", "31", "35", "37", "41", "54", "84", "100", "112", "116", "118", "133", "154", "199", "203", "204", "208", "219", "294", "318", "339", "353", "363", "378")
adult.treated.high=c("21", "22", "45", "82", "85", "91", "95", "99", "103", "108", "117", "127", "165", "185", "190", "196", "232", "237", "245", "263", "276", "306", "343", "374")
# Combine vectors into lists
# Used for adding treatment info and/or subsetting downstream
groups_list <- list(juvenile.control.ambient = juvenile.control.ambient,
juvenile.control.high = juvenile.control.high,
juvenile.treated.ambient = juvenile.treated.ambient,
juvenile.treated.high = juvenile.treated.high,
adult.control.ambient = adult.control.ambient,
adult.control.high = adult.control.high,
adult.treated.ambient = adult.treated.ambient,
adult.treated.high = adult.treated.high,
seed.control.ambient = seed.control.ambient,
seed.control.high = seed.control.high,
seed.treated.ambient = seed.treated.ambient,
seed.treated.high = seed.treated.high,
spat.control.ambient = spat.control.ambient,
spat.control.high = spat.control.high,
spat.treated.ambient = spat.treated.ambient,
spat.treated.high = spat.treated.high)
Normalized to designated normalizing gene
calculate_delta_Cq <- function(df) {
df <- df %>%
group_by(Sample) %>%
mutate(delta_Cq = Cq.Mean - Cq.Mean[Target == "GAPDH"]) %>%
ungroup()
return(df)
}
# Function to create box plots for each comparison
create_boxplot_delta_Cq <- function(data, comparison, t_test_results) {
# Extract life stages from comparison
life_stages <- unlist(strsplit(comparison, "\\."))
# Debugging: Print life stages
# print(paste("Life stages for comparison:", comparison))
# print(life_stages)
# Filter data for the relevant life stages
filtered_data <- data %>%
filter(life.stage %in% life_stages)
# Debugging: Print filtered data
# print("Filtered data:")
# print(filtered_data)
# Check if both life stages are included
if (!all(life_stages %in% unique(filtered_data$life.stage))) {
stop("Not all life stages are included in the filtered data")
}
y_limits <- range(filtered_data$delta_Cq, na.rm = TRUE)
# Debugging: Print y_limits
# print("Y limits:")
# print(y_limits)
# Filter t_test_results for the current comparison
t_test_results_filtered <- t_test_results %>%
filter(comparison == !!comparison)
# Debugging: Print filtered t_test_results
# print("Filtered t_test_results:")
# print(t_test_results_filtered)
# Filter t_test_results for asterisks
t_test_results_with_asterisks <- t_test_results_filtered %>%
filter(asterisk != "")
# Debugging: Print t_test_results_with_asterisks
# print("t_test_results_with_asterisks:")
# print(t_test_results_with_asterisks)
formatted_title <- paste0(toupper(substring(life_stages[1], 1, 1)), substring(life_stages[1], 2),
" vs. ",
toupper(substring(life_stages[2], 1, 1)), substring(life_stages[2], 2))
boxplot <- ggplot(filtered_data, aes(x = Target, y = delta_Cq, fill = life.stage)) +
geom_boxplot(position = position_dodge(width = 0.75)) +
theme_minimal() +
theme(legend.position = "right") +
scale_fill_manual(values=c("darkgray", "salmon", "lightblue", "lightgreen")) +
ylim(y_limits) +
labs(x = "Target", y = "Delta Cq", title = formatted_title) +
# Highlighted section: Adds asterisks
geom_text(data = t_test_results_with_asterisks,
aes(x = Target, y = y_limits[2] - 1, label = asterisk),
vjust = -0.5, size = 8, color = "black", inherit.aes = FALSE)
print(boxplot)
}
life_stage, treatment1, and treatment2).filtered_data.The t_test_results_filtered data frame is filtered for the specific comparison.
The t_test_results_with_asterisks data frame is created to include only the rows with asterisks.
The formatted_title variable is created by capitalizing the first letter of each component and concatenating them with ” - ” and ” vs. ” in between.
This should create box plots comparing conditioning treatments within each life stage, with titles formatted as <life.stage> - Treated vs. Control.
# Function to create box plots for each comparison of conditioning treatments within life stages
create_boxplot_conditioning <- function(data, comparison, t_test_results) {
# Extract life stage and conditioning treatments from comparison
comparison_parts <- unlist(strsplit(comparison, "\\."))
life_stage <- comparison_parts[1]
treatment1 <- comparison_parts[2]
treatment2 <- comparison_parts[3]
# Debugging: Print life stage and treatments
# print(paste("Life stage and treatments for comparison:", comparison))
# print(c(life_stage, treatment1, treatment2))
# Filter data for the relevant life stage and conditioning treatments
filtered_data <- data %>%
filter(life.stage == life_stage, conditioning.treatment %in% c(treatment1, treatment2))
# Debugging: Print filtered data
# print("Filtered data:")
# print(filtered_data)
# Check if both treatments are included
if (!all(c(treatment1, treatment2) %in% unique(filtered_data$conditioning.treatment))) {
stop("Not all treatments are included in the filtered data")
}
y_limits <- range(filtered_data$delta_Cq, na.rm = TRUE)
# Debugging: Print y_limits
# print("Y limits:")
# print(y_limits)
# Filter t_test_results for the current comparison
t_test_results_filtered <- t_test_results %>%
filter(comparison == !!comparison)
# Debugging: Print filtered t_test_results
# print("Filtered t_test_results:")
# print(t_test_results_filtered)
# Filter t_test_results for asterisks
t_test_results_with_asterisks <- t_test_results_filtered %>%
filter(asterisk != "")
# Debugging: Print t_test_results_with_asterisks
# print("t_test_results_with_asterisks:")
# print(t_test_results_with_asterisks)
# Format the title
formatted_title <- paste0(toupper(substring(life_stage, 1, 1)), substring(life_stage, 2),
" - ",
toupper(substring(treatment1, 1, 1)), substring(treatment1, 2),
" vs. ",
toupper(substring(treatment2, 1, 1)), substring(treatment2, 2))
boxplot <- ggplot(filtered_data, aes(x = Target, y = delta_Cq, fill = conditioning.treatment)) +
geom_boxplot(position = position_dodge(width = 0.75)) +
theme_minimal() +
theme(legend.position = "right") +
scale_fill_manual(values=c("darkgray", "salmon")) +
ylim(y_limits) +
labs(x = "Target", y = "Delta Cq", title = formatted_title) +
# Highlighted section: Adds asterisks
geom_text(data = t_test_results_with_asterisks,
aes(x = Target, y = y_limits[2] - 1, label = asterisk),
vjust = -0.5, size = 8, color = "black", inherit.aes = FALSE)
print(boxplot)
}
life_stage, treatment1, and treatment2).filtered_data data frame is filtered to include only the rows with the relevant life stage and acute treatments.filtered_data.The t_test_results_filtered data frame is filtered for the specific comparison.
The t_test_results_with_asterisks data frame is created to include only the rows with asterisks.
Format the Title:
<life.stage> - Ambient vs. High.# Function to create box plots for each comparison of acute treatments within life stages
create_boxplot_acute <- function(data, comparison, t_test_results) {
# Extract life stage and acute treatments from comparison
comparison_parts <- unlist(strsplit(comparison, "\\."))
life_stage <- comparison_parts[1]
treatment1 <- comparison_parts[2]
treatment2 <- comparison_parts[3]
# Debugging: Print life stage and treatments
# print(paste("Life stage and treatments for comparison:", comparison))
# print(c(life_stage, treatment1, treatment2))
# Filter data for the relevant life stage and acute treatments
filtered_data <- data %>%
filter(life.stage == life_stage, acute.treatment %in% c(treatment1, treatment2))
# Debugging: Print filtered data
# print("Filtered data:")
# print(filtered_data)
# Check if both treatments are included
if (!all(c(treatment1, treatment2) %in% unique(filtered_data$acute.treatment))) {
stop("Not all treatments are included in the filtered data")
}
y_limits <- range(filtered_data$delta_Cq, na.rm = TRUE)
# Debugging: Print y_limits
# print("Y limits:")
# print(y_limits)
# Filter t_test_results for the current comparison
t_test_results_filtered <- t_test_results %>%
filter(comparison == !!comparison)
# Debugging: Print filtered t_test_results
# print("Filtered t_test_results:")
# print(t_test_results_filtered)
# Filter t_test_results for asterisks
t_test_results_with_asterisks <- t_test_results_filtered %>%
filter(asterisk != "")
# Debugging: Print t_test_results_with_asterisks
# print("t_test_results_with_asterisks:")
# print(t_test_results_with_asterisks)
# Format the title
formatted_title <- paste0(toupper(substring(life_stage, 1, 1)), substring(life_stage, 2),
" - ",
toupper(substring(treatment1, 1, 1)), substring(treatment1, 2),
" vs. ",
toupper(substring(treatment2, 1, 1)), substring(treatment2, 2))
boxplot <- ggplot(filtered_data, aes(x = Target, y = delta_Cq, fill = acute.treatment)) +
geom_boxplot(position = position_dodge(width = 0.75)) +
theme_minimal() +
theme(legend.position = "right") +
scale_fill_manual(values=c("darkgray", "salmon")) +
ylim(y_limits) +
labs(x = "Target", y = "Delta Cq", title = formatted_title) +
# Highlighted section: Adds asterisks
geom_text(data = t_test_results_with_asterisks,
aes(x = Target, y = y_limits[2] - 1, label = asterisk),
vjust = -0.5, size = 8, color = "magenta", inherit.aes = FALSE)
print(boxplot)
}
# Function to create box plots for each comparison of acute treatments within life stages and conditioning treatments
create_boxplot_acute_conditioning <- function(data, comparison, t_test_results) {
# Extract life stage, conditioning treatment, and acute treatments from comparison
comparison_parts <- unlist(strsplit(comparison, "\\."))
life_stage <- comparison_parts[1]
conditioning_treatment <- comparison_parts[2]
treatment1 <- comparison_parts[3]
treatment2 <- comparison_parts[5]
# Filter data for the relevant life stage, conditioning treatment, and acute treatments
filtered_data <- data %>%
filter(life.stage == life_stage, conditioning.treatment == conditioning_treatment, acute.treatment %in% c(treatment1, treatment2))
# Check if both treatments are included
if (!all(c(treatment1, treatment2) %in% unique(filtered_data$acute.treatment))) {
stop("Not all treatments are included in the filtered data")
}
y_limits <- range(filtered_data$delta_Cq, na.rm = TRUE)
# Filter t_test_results for the current comparison
t_test_results_filtered <- t_test_results %>%
filter(comparison == !!comparison)
# Filter t_test_results for asterisks
t_test_results_with_asterisks <- t_test_results_filtered %>%
filter(asterisk != "")
# Format the title
formatted_title <- paste0(toupper(substring(life_stage, 1, 1)), substring(life_stage, 2),
" - ",
toupper(substring(conditioning_treatment, 1, 1)), substring(conditioning_treatment, 2),
" - ",
toupper(substring(treatment1, 1, 1)), substring(treatment1, 2),
" vs. ",
toupper(substring(treatment2, 1, 1)), substring(treatment2, 2))
boxplot <- ggplot(filtered_data, aes(x = Target, y = delta_Cq, fill = acute.treatment)) +
geom_boxplot(position = position_dodge(width = 0.75)) +
theme_minimal() +
theme(legend.position = "right") +
scale_fill_manual(values=c("darkgray", "salmon")) +
ylim(y_limits) +
labs(x = "Target", y = "Delta Cq", title = formatted_title) +
# Adds asterisks
geom_text(data = t_test_results_with_asterisks,
aes(x = Target, y = y_limits[2] - 1, label = asterisk),
vjust = -0.5, size = 8, color = "black", inherit.aes = FALSE)
print(boxplot)
}
# Get a list of all CSV files in the directory with the naming structure "*Cq-Results.csv"
cq_file_list <- list() # Initialize list
cq_file_list <- list.files(path = cqs_directory, pattern = "Cq-Results\\.csv$", full.names = TRUE)
# Initialize an empty list to store the data frames
data_frames_list <- list()
# Loop through each file and read it into a data frame, then add it to the list
for (file in cq_file_list) {
data <- read.csv(file, header = TRUE)
data$Sample <- as.character(data$Sample) # Convert Sample column to character type
data_frames_list[[file]] <- data
}
# Combine all data frames into a single data frame
combined_df <- bind_rows(data_frames_list, .id = "data_frame_id")
str(combined_df)
'data.frame': 2192 obs. of 17 variables:
$ data_frame_id : chr "lifestage_carryover/data/qPCR/Cq/sam_2024-03-25_06-10-54_Connect-Quantification-Cq-Results.csv" "lifestage_carryover/data/qPCR/Cq/sam_2024-03-25_06-10-54_Connect-Quantification-Cq-Results.csv" "lifestage_carryover/data/qPCR/Cq/sam_2024-03-25_06-10-54_Connect-Quantification-Cq-Results.csv" "lifestage_carryover/data/qPCR/Cq/sam_2024-03-25_06-10-54_Connect-Quantification-Cq-Results.csv" ...
$ X : logi NA NA NA NA NA NA ...
$ Well : chr "A01" "A02" "A03" "A04" ...
$ Fluor : chr "SYBR" "SYBR" "SYBR" "SYBR" ...
$ Target : chr "ATPsynthase" "ATPsynthase" "ATPsynthase" "ATPsynthase" ...
$ Content : chr "Unkn-01" "Unkn-01" "Unkn-02" "Unkn-02" ...
$ Sample : chr "206" "206" "220" "220" ...
$ Biological.Set.Name : logi NA NA NA NA NA NA ...
$ Cq : num 26.7 26.7 25.8 25.9 25.1 ...
$ Cq.Mean : num 26.7 26.7 25.9 25.9 25.1 ...
$ Cq.Std..Dev : num 0.0455 0.0455 0.0239 0.0239 0.0813 ...
$ Starting.Quantity..SQ.: num NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
$ Log.Starting.Quantity : num NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
$ SQ.Mean : num NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
$ SQ.Std..Dev : num NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
$ Set.Point : int 60 60 60 60 60 60 60 60 60 60 ...
$ Well.Note : logi NA NA NA NA NA NA ...
# Remove rows with Sample name "NTC"
combined_df <- combined_df[combined_df$Sample != "NTC", ]
# Replace values in the Target column
combined_df$Target <- gsub("Cg_GAPDH_205_F-355_R \\(SR IDs: 1172/3\\)", "GAPDH", combined_df$Target)
combined_df$Target <- gsub("Cg_ATPsynthase_F/R \\(SR IDs: 1385/6\\)", "ATPsynthase", combined_df$Target)
combined_df$Target <- gsub("Cg_cGAS \\(SR IDs: 1826/7\\)", "cGAS", combined_df$Target)
combined_df$Target <- gsub("Cg_citrate_synthase \\(SR IDs: 1383/4\\)", "citrate synthase", combined_df$Target)
combined_df$Target <- gsub("Cg_DNMT1_F \\(SR IDs: 1510/1\\)", "DNMT1", combined_df$Target)
combined_df$Target <- gsub("Cg_HSP70_F/R \\(SR IDs: 598/9\\)", "HSP70", combined_df$Target)
combined_df$Target <- gsub("Cg_Hsp90_F/R \\(SR IDs: 1532/3\\)", "HSP90", combined_df$Target)
combined_df$Target <- gsub("Cg_VIPERIN_F/R \\(SR IDs: 1828/9\\)", "viperin", combined_df$Target)
str(combined_df)
'data.frame': 2180 obs. of 17 variables:
$ data_frame_id : chr "lifestage_carryover/data/qPCR/Cq/sam_2024-03-25_06-10-54_Connect-Quantification-Cq-Results.csv" "lifestage_carryover/data/qPCR/Cq/sam_2024-03-25_06-10-54_Connect-Quantification-Cq-Results.csv" "lifestage_carryover/data/qPCR/Cq/sam_2024-03-25_06-10-54_Connect-Quantification-Cq-Results.csv" "lifestage_carryover/data/qPCR/Cq/sam_2024-03-25_06-10-54_Connect-Quantification-Cq-Results.csv" ...
$ X : logi NA NA NA NA NA NA ...
$ Well : chr "A01" "A02" "A03" "A04" ...
$ Fluor : chr "SYBR" "SYBR" "SYBR" "SYBR" ...
$ Target : chr "ATPsynthase" "ATPsynthase" "ATPsynthase" "ATPsynthase" ...
$ Content : chr "Unkn-01" "Unkn-01" "Unkn-02" "Unkn-02" ...
$ Sample : chr "206" "206" "220" "220" ...
$ Biological.Set.Name : logi NA NA NA NA NA NA ...
$ Cq : num 26.7 26.7 25.8 25.9 25.1 ...
$ Cq.Mean : num 26.7 26.7 25.9 25.9 25.1 ...
$ Cq.Std..Dev : num 0.0455 0.0455 0.0239 0.0239 0.0813 ...
$ Starting.Quantity..SQ.: num NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
$ Log.Starting.Quantity : num NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
$ SQ.Mean : num NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
$ SQ.Std..Dev : num NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
$ Set.Point : int 60 60 60 60 60 60 60 60 60 60 ...
$ Well.Note : logi NA NA NA NA NA NA ...
levels(as.factor(combined_df$Target))
[1] "ATPsynthase" "cGAS" "citrate synthase" "DNMT1"
[5] "GAPDH" "HSP70" "HSP90" "viperin"
# Filter out rows where Cq.Std..Dev is NA
combined_df <- combined_df[!is.na(combined_df$Cq.Std..Dev), ]
# Filter rows where Cq.Std..Dev is greater than 0.5
high_cq_std_dev <- combined_df[combined_df$Cq.Std..Dev > 0.5, ]
# Print the filtered rows with specified columns, without row names
print(high_cq_std_dev[, c("Target", "Sample", "Cq", "Cq.Std..Dev")], row.names = FALSE)
Target Sample Cq Cq.Std..Dev
HSP70 244 33.10339 4.0838809
HSP70 244 27.32791 4.0838809
HSP90 223 24.85319 0.7548714
HSP90 223 25.92074 0.7548714
viperin 223 30.30089 0.6058663
viperin 223 31.15772 0.6058663
viperin 243 32.57817 0.5527617
viperin 243 33.35989 0.5527617
DNMT1 296 31.21374 0.6417578
DNMT1 296 30.30616 0.6417578
DNMT1 298 35.68716 0.5406704
DNMT1 298 34.92253 0.5406704
DNMT1 223 32.24089 0.6214201
DNMT1 223 33.11971 0.6214201
DNMT1 243 36.63921 0.5125743
DNMT1 243 35.91432 0.5125743
DNMT1 285 33.63443 0.7036122
DNMT1 285 34.62949 0.7036122
GAPDH 316 23.94926 8.5684728
GAPDH 316 24.14183 8.5684728
GAPDH 316 38.88564 8.5684728
GAPDH 213 26.98012 2.2910353
GAPDH 213 23.00009 2.2910353
GAPDH 213 26.95634 2.2910353
GAPDH 263 22.42154 0.8731474
GAPDH 263 23.77008 0.8731474
GAPDH 263 24.05667 0.8731474
citrate synthase 230 24.44066 4.4783429
citrate synthase 230 24.40421 4.4783429
citrate synthase 230 32.17909 4.4783429
viperin 227 30.47773 3.5152533
viperin 227 30.37738 3.5152533
viperin 227 36.51553 3.5152533
viperin 245 26.05748 5.1635899
viperin 245 34.98192 5.1635899
viperin 245 26.01928 5.1635899
viperin 341 26.48675 2.9838590
viperin 341 31.67235 2.9838590
viperin 341 26.52174 2.9838590
viperin 344 29.98184 2.3712440
viperin 344 25.90358 2.3712440
viperin 344 25.84648 2.3712440
viperin 355 28.79712 0.5821437
viperin 355 29.57428 0.5821437
viperin 355 28.43490 0.5821437
# Group by Sample and Target, then filter out the outlier replicate
combined.fitered_df<- combined_df %>%
group_by(Sample, Target) %>%
filter(abs(Cq - mean(Cq, na.rm = TRUE)) <= Cq.Std..Dev)
# Print the filtered data frame
str(combined.fitered_df)
gropd_df [1,520 × 17] (S3: grouped_df/tbl_df/tbl/data.frame)
$ data_frame_id : chr [1:1520] "lifestage_carryover/data/qPCR/Cq/sam_2024-03-25_06-10-54_Connect-Quantification-Cq-Results.csv" "lifestage_carryover/data/qPCR/Cq/sam_2024-03-25_06-10-54_Connect-Quantification-Cq-Results.csv" "lifestage_carryover/data/qPCR/Cq/sam_2024-03-25_06-10-54_Connect-Quantification-Cq-Results.csv" "lifestage_carryover/data/qPCR/Cq/sam_2024-03-25_06-10-54_Connect-Quantification-Cq-Results.csv" ...
$ X : logi [1:1520] NA NA NA NA NA NA ...
$ Well : chr [1:1520] "A01" "A02" "A03" "A04" ...
$ Fluor : chr [1:1520] "SYBR" "SYBR" "SYBR" "SYBR" ...
$ Target : chr [1:1520] "ATPsynthase" "ATPsynthase" "ATPsynthase" "ATPsynthase" ...
$ Content : chr [1:1520] "Unkn-01" "Unkn-01" "Unkn-02" "Unkn-02" ...
$ Sample : chr [1:1520] "206" "206" "220" "220" ...
$ Biological.Set.Name : logi [1:1520] NA NA NA NA NA NA ...
$ Cq : num [1:1520] 26.7 26.7 25.8 25.9 25.1 ...
$ Cq.Mean : num [1:1520] 26.7 26.7 25.9 25.9 25.1 ...
$ Cq.Std..Dev : num [1:1520] 0.0455 0.0455 0.0239 0.0239 0.0813 ...
$ Starting.Quantity..SQ.: num [1:1520] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
$ Log.Starting.Quantity : num [1:1520] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
$ SQ.Mean : num [1:1520] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
$ SQ.Std..Dev : num [1:1520] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
$ Set.Point : int [1:1520] 60 60 60 60 60 60 60 60 60 60 ...
$ Well.Note : logi [1:1520] NA NA NA NA NA NA ...
- attr(*, "groups")= tibble [760 × 3] (S3: tbl_df/tbl/data.frame)
..$ Sample: chr [1:760] "200" "200" "200" "200" ...
..$ Target: chr [1:760] "ATPsynthase" "DNMT1" "GAPDH" "HSP70" ...
..$ .rows : list<int> [1:760]
.. ..$ : int [1:2] 21 22
.. ..$ : int [1:2] 245 246
.. ..$ : int [1:2] 85 86
.. ..$ : int [1:2] 53 54
.. ..$ : int [1:2] 117 118
.. ..$ : int [1:2] 149 150
.. ..$ : int [1:2] 213 214
.. ..$ : int [1:2] 181 182
.. ..$ : int [1:2] 321 322
.. ..$ : int [1:2] 889 890
.. ..$ : int [1:2] 637 638
.. ..$ : int [1:2] 1047 1048
.. ..$ : int [1:2] 1205 1206
.. ..$ : int [1:2] 479 480
.. ..$ : int [1:2] 731 732
.. ..$ : int [1:2] 1363 1364
.. ..$ : int [1:2] 323 324
.. ..$ : int [1:2] 891 892
.. ..$ : int [1:2] 639 640
.. ..$ : int [1:2] 1049 1050
.. ..$ : int [1:2] 1207 1208
.. ..$ : int [1:2] 481 482
.. ..$ : int [1:2] 733 734
.. ..$ : int [1:2] 1365 1366
.. ..$ : int [1:2] 325 326
.. ..$ : int [1:2] 893 894
.. ..$ : int [1:2] 641 642
.. ..$ : int [1:2] 1051 1052
.. ..$ : int [1:2] 1209 1210
.. ..$ : int [1:2] 483 484
.. ..$ : int [1:2] 735 736
.. ..$ : int [1:2] 1367 1368
.. ..$ : int [1:2] 327 328
.. ..$ : int [1:2] 895 896
.. ..$ : int [1:2] 643 644
.. ..$ : int [1:2] 1053 1054
.. ..$ : int [1:2] 1211 1212
.. ..$ : int [1:2] 485 486
.. ..$ : int [1:2] 737 738
.. ..$ : int [1:2] 1369 1370
.. ..$ : int [1:2] 329 330
.. ..$ : int [1:2] 897 898
.. ..$ : int [1:2] 645 646
.. ..$ : int [1:2] 1055 1056
.. ..$ : int [1:2] 1213 1214
.. ..$ : int [1:2] 487 488
.. ..$ : int [1:2] 739 740
.. ..$ : int [1:2] 1371 1372
.. ..$ : int [1:2] 1 2
.. ..$ : int [1:2] 225 226
.. ..$ : int [1:2] 65 66
.. ..$ : int [1:2] 33 34
.. ..$ : int [1:2] 97 98
.. ..$ : int [1:2] 129 130
.. ..$ : int [1:2] 193 194
.. ..$ : int [1:2] 161 162
.. ..$ : int [1:2] 331 332
.. ..$ : int [1:2] 899 900
.. ..$ : int [1:2] 647 648
.. ..$ : int [1:2] 1057 1058
.. ..$ : int [1:2] 1215 1216
.. ..$ : int [1:2] 489 490
.. ..$ : int [1:2] 741 742
.. ..$ : int [1:2] 1373 1374
.. ..$ : int [1:2] 333 334
.. ..$ : int [1:2] 901 902
.. ..$ : int [1:2] 649 650
.. ..$ : int [1:2] 1059 1060
.. ..$ : int [1:2] 1217 1218
.. ..$ : int [1:2] 491 492
.. ..$ : int [1:2] 743 744
.. ..$ : int [1:2] 1375 1376
.. ..$ : int [1:2] 335 336
.. ..$ : int [1:2] 903 904
.. ..$ : int [1:2] 651 652
.. ..$ : int [1:2] 1061 1062
.. ..$ : int [1:2] 1219 1220
.. ..$ : int [1:2] 493 494
.. ..$ : int [1:2] 745 746
.. ..$ : int [1:2] 1377 1378
.. ..$ : int [1:2] 337 338
.. ..$ : int [1:2] 905 906
.. ..$ : int [1:2] 653 654
.. ..$ : int [1:2] 1063 1064
.. ..$ : int [1:2] 1221 1222
.. ..$ : int [1:2] 495 496
.. ..$ : int [1:2] 747 748
.. ..$ : int [1:2] 1379 1380
.. ..$ : int [1:2] 339 340
.. ..$ : int [1:2] 907 908
.. ..$ : int [1:2] 655 656
.. ..$ : int [1:2] 1065 1066
.. ..$ : int [1:2] 1223 1224
.. ..$ : int [1:2] 497 498
.. ..$ : int [1:2] 749 750
.. ..$ : int [1:2] 1381 1382
.. ..$ : int [1:2] 341 342
.. ..$ : int [1:2] 909 910
.. ..$ : int [1:2] 657 658
.. .. [list output truncated]
.. ..@ ptype: int(0)
..- attr(*, ".drop")= logi TRUE
# Group by Sample and Target, then summarize to get unique rows for each sample
grouped_df <- combined.fitered_df%>%
group_by(Sample, Target) %>%
summarize(Cq.Mean = mean(Cq, na.rm = TRUE)) %>%
ungroup()
str(grouped_df)
tibble [760 × 3] (S3: tbl_df/tbl/data.frame)
$ Sample : chr [1:760] "200" "200" "200" "200" ...
$ Target : chr [1:760] "ATPsynthase" "DNMT1" "GAPDH" "HSP70" ...
$ Cq.Mean: num [1:760] 25.2 33.6 25.4 31.8 26 ...
# Initialize new columns
grouped_df <- grouped_df %>%
mutate(life.stage = NA_character_,
conditioning.treatment = NA_character_,
acute.treatment = NA_character_)
# Loop through each vector
for (vec_name in names(groups_list)) {
vec <- groups_list[[vec_name]]
stage <- strsplit(vec_name, "\\.")[[1]][1]
conditioning_treatment <- strsplit(vec_name, "\\.")[[1]][2]
acute_treatment <- strsplit(vec_name, "\\.")[[1]][3]
# Loop through each row in grouped_df
for (i in 1:nrow(grouped_df)) {
sample <- grouped_df$Sample[i]
# Check if sample is in the vector
if (sample %in% vec) {
# Update life.stage and treatment columns
grouped_df$life.stage[i] <- stage
grouped_df$conditioning.treatment[i] <- conditioning_treatment
grouped_df$acute.treatment[i] <-acute_treatment
}
}
}
str(grouped_df)
tibble [760 × 6] (S3: tbl_df/tbl/data.frame)
$ Sample : chr [1:760] "200" "200" "200" "200" ...
$ Target : chr [1:760] "ATPsynthase" "DNMT1" "GAPDH" "HSP70" ...
$ Cq.Mean : num [1:760] 25.2 33.6 25.4 31.8 26 ...
$ life.stage : chr [1:760] "seed" "seed" "seed" "seed" ...
$ conditioning.treatment: chr [1:760] "treated" "treated" "treated" "treated" ...
$ acute.treatment : chr [1:760] "high" "high" "high" "high" ...
# Calculate delta Cq by subtracting GAPDH Cq.Mean from each corresponding Sample Cq.Mean
delta_Cq_df <- calculate_delta_Cq(grouped_df)
# Filters out normalizing gene, since no need to compare normalizing gene to itself.
delta_Cq_df <- delta_Cq_df %>%
filter(!is.na(life.stage), !is.na(Target), Target != "GAPDH")
str(delta_Cq_df)
tibble [665 × 7] (S3: tbl_df/tbl/data.frame)
$ Sample : chr [1:665] "200" "200" "200" "200" ...
$ Target : chr [1:665] "ATPsynthase" "DNMT1" "HSP70" "HSP90" ...
$ Cq.Mean : num [1:665] 25.2 33.6 31.8 26 30.6 ...
$ life.stage : chr [1:665] "seed" "seed" "seed" "seed" ...
$ conditioning.treatment: chr [1:665] "treated" "treated" "treated" "treated" ...
$ acute.treatment : chr [1:665] "high" "high" "high" "high" ...
$ delta_Cq : num [1:665] -0.243 8.168 6.349 0.578 5.135 ...
This code does the following:
# Extract unique life.stage levels
unique_life_stages <- unique(delta_Cq_df$life.stage)
# Generate all possible pairs of life.stage levels
life_stage_pairs <- combn(unique_life_stages, 2, simplify = FALSE)
# Initialize a list to store results
life_stage_t_test_results_list <- list()
for (pair in life_stage_pairs) {
stage1 <- pair[1]
stage2 <- pair[2]
# Perform t-test for each Target comparing the two life.stage levels
t_test_results <- delta_Cq_df %>%
filter(life.stage %in% c(stage1, stage2)) %>%
group_by(Target) %>%
summarise(
t_test_result = list(t.test(delta_Cq ~ life.stage))
) %>%
ungroup() %>%
mutate(
estimate_diff = sapply(t_test_result, function(x) x$estimate[1] - x$estimate[2]),
p_value = sapply(t_test_result, function(x) x$p.value),
asterisk = ifelse(p_value <= 0.05, "*", ""), # Adds asterisk column and asterisk for p-value.
comparison = paste(stage1, "vs", stage2, sep = ".")
) %>%
select(!t_test_result)
life_stage_t_test_results_list[[paste(stage1, stage2, sep = ".")]] <- t_test_results
}
# Combine results into a single data frame
life_stage_t_test_results_df <- bind_rows(life_stage_t_test_results_list, .id = "comparison")
# View the results
print(life_stage_t_test_results_df)
# A tibble: 42 × 5
Target estimate_diff p_value asterisk comparison
<chr> <dbl> <dbl> <chr> <chr>
1 ATPsynthase 0.310 0.0183 "*" seed.adult
2 DNMT1 -0.0736 0.829 "" seed.adult
3 HSP70 0.179 0.702 "" seed.adult
4 HSP90 0.727 0.00635 "*" seed.adult
5 cGAS -0.0136 0.963 "" seed.adult
6 citrate synthase -0.215 0.391 "" seed.adult
7 viperin -1.08 0.0000737 "*" seed.adult
8 ATPsynthase 0.207 0.0819 "" seed.juvenile
9 DNMT1 -0.434 0.254 "" seed.juvenile
10 HSP70 0.657 0.179 "" seed.juvenile
# ℹ 32 more rows
This code does the following:
# Extract unique life.stage levels
unique_life_stages <- unique(delta_Cq_df$life.stage)
# Initialize a list to store results
conditioning_treatment_t_test_results_list <- list()
for (stage in unique_life_stages) {
# Extract unique conditioning.treatment levels within the current life.stage
unique_treatments <- unique(delta_Cq_df %>% filter(life.stage == stage) %>% pull(conditioning.treatment))
# Generate all possible pairs of conditioning.treatment levels
treatment_pairs <- combn(unique_treatments, 2, simplify = FALSE)
for (pair in treatment_pairs) {
treatment1 <- pair[1]
treatment2 <- pair[2]
# Perform t-test for each Target comparing the two conditioning.treatment levels within the current life.stage
t_test_results <- delta_Cq_df %>%
filter(life.stage == stage, conditioning.treatment %in% c(treatment1, treatment2)) %>%
group_by(Target) %>%
summarise(
t_test_result = list(t.test(delta_Cq ~ conditioning.treatment))
) %>%
ungroup() %>%
mutate(
estimate_diff = sapply(t_test_result, function(x) x$estimate[1] - x$estimate[2]),
p_value = sapply(t_test_result, function(x) x$p.value),
asterisk = ifelse(p_value <= 0.05, "*", ""), # Adds asterisk column and asterisk for p-value.
comparison = paste(stage, treatment1, "vs", treatment2, sep = ".")
) %>%
select(!t_test_result)
conditioning_treatment_t_test_results_list[[paste(stage, treatment1, treatment2, sep = ".")]] <- t_test_results
}
}
# Combine results into a single data frame
conditioning_treatment_t_test_results_df <- bind_rows(conditioning_treatment_t_test_results_list, .id = "comparison")
# View the results
print(conditioning_treatment_t_test_results_df)
# A tibble: 28 × 5
Target estimate_diff p_value asterisk comparison
<chr> <dbl> <dbl> <chr> <chr>
1 ATPsynthase -0.0248 0.907 "" seed.treated.control
2 DNMT1 0.206 0.746 "" seed.treated.control
3 HSP70 -0.123 0.851 "" seed.treated.control
4 HSP90 -0.299 0.439 "" seed.treated.control
5 cGAS -0.0776 0.888 "" seed.treated.control
6 citrate synthase -0.0148 0.976 "" seed.treated.control
7 viperin 0.160 0.695 "" seed.treated.control
8 ATPsynthase 0.0779 0.603 "" adult.treated.control
9 DNMT1 0.312 0.278 "" adult.treated.control
10 HSP70 -0.941 0.177 "" adult.treated.control
# ℹ 18 more rows
This code does the following:
Excludes seed and spat, as these were only held at ambient for the acute treatment.
# Extract unique life.stage levels, excluding 'seed' and 'spat'
unique_life_stages <- unique(delta_Cq_df$life.stage)
unique_life_stages <- setdiff(unique_life_stages, c("seed", "spat"))
# Initialize a list to store results
acute_treatment_t_test_results_list <- list()
for (stage in unique_life_stages) {
# Extract unique acute.treatment levels within the current life.stage
unique_treatments <- unique(delta_Cq_df %>% filter(life.stage == stage) %>% pull(acute.treatment))
# Check if there are at least 2 unique treatments
if (length(unique_treatments) >= 2) {
# Generate all possible pairs of acute.treatment levels
treatment_pairs <- combn(unique_treatments, 2, simplify = FALSE)
for (pair in treatment_pairs) {
treatment1 <- pair[1]
treatment2 <- pair[2]
# Perform t-test for each Target comparing the two acute.treatment levels within the current life.stage
t_test_results <- delta_Cq_df %>%
filter(life.stage == stage, acute.treatment %in% c(treatment1, treatment2)) %>%
group_by(Target) %>%
summarise(
t_test_result = list(t.test(delta_Cq ~ acute.treatment))
) %>%
ungroup() %>%
mutate(
estimate_diff = sapply(t_test_result, function(x) x$estimate[1] - x$estimate[2]),
p_value = sapply(t_test_result, function(x) x$p.value),
asterisk = ifelse(p_value <= 0.05, "*", ""), # Adds asterisk column and asterisk for p-value.
comparison = paste(stage, treatment1, "vs", treatment2, sep = ".")
) %>%
select(!t_test_result)
acute_treatment_t_test_results_list[[paste(stage, treatment1, treatment2, sep = ".")]] <- t_test_results
}
}
}
# Combine results into a single data frame
acute_treatment_t_test_results_df <- bind_rows(acute_treatment_t_test_results_list, .id = "comparison")
# View the results
print(acute_treatment_t_test_results_df)
# A tibble: 14 × 5
Target estimate_diff p_value asterisk comparison
<chr> <dbl> <dbl> <chr> <chr>
1 ATPsynthase 0.0605 0.687 "" adult.ambient.high
2 DNMT1 0.314 0.275 "" adult.ambient.high
3 HSP70 0.276 0.696 "" adult.ambient.high
4 HSP90 0.499 0.149 "" adult.ambient.high
5 cGAS 0.329 0.200 "" adult.ambient.high
6 citrate synthase 0.0668 0.622 "" adult.ambient.high
7 viperin 0.323 0.251 "" adult.ambient.high
8 ATPsynthase -0.0370 0.738 "" juvenile.ambient.high
9 DNMT1 -0.672 0.121 "" juvenile.ambient.high
10 HSP70 0.745 0.319 "" juvenile.ambient.high
11 HSP90 0.0450 0.859 "" juvenile.ambient.high
12 cGAS -0.203 0.474 "" juvenile.ambient.high
13 citrate synthase 0.0399 0.870 "" juvenile.ambient.high
14 viperin -0.424 0.304 "" juvenile.ambient.high
# Extract unique life.stage levels, excluding 'seed' and 'spat'
unique_life_stages <- unique(delta_Cq_df$life.stage)
#unique_life_stages <- setdiff(unique_life_stages, c("seed", "spat"))
# Extract unique conditioning.treatment levels
unique_conditioning_treatments <- unique(delta_Cq_df$conditioning.treatment)
# Initialize a list to store results
acute_treatment_within_life.stages_conditioning_t_test_results_list <- list()
for (stage in unique_life_stages) {
for (conditioning in unique_conditioning_treatments) {
# Extract unique acute.treatment levels within the current life.stage and conditioning.treatment
unique_treatments <- unique(delta_Cq_df %>% filter(life.stage == stage, conditioning.treatment == conditioning) %>% pull(acute.treatment))
# Check if there are at least 2 unique treatments
if (length(unique_treatments) >= 2) {
# Generate all possible pairs of acute.treatment levels
treatment_pairs <- combn(unique_treatments, 2, simplify = FALSE)
for (pair in treatment_pairs) {
treatment1 <- pair[1]
treatment2 <- pair[2]
# Perform t-test for each Target comparing the two acute.treatment levels within the current life.stage and conditioning.treatment
t_test_results <- delta_Cq_df %>%
filter(life.stage == stage, conditioning.treatment == conditioning, acute.treatment %in% c(treatment1, treatment2)) %>%
group_by(Target) %>%
summarise(
t_test_result = list(t.test(delta_Cq ~ acute.treatment))
) %>%
ungroup() %>%
mutate(
estimate_diff = sapply(t_test_result, function(x) x$estimate[1] - x$estimate[2]),
p_value = sapply(t_test_result, function(x) x$p.value),
asterisk = ifelse(p_value <= 0.05, "*", ""), # Adds asterisk column and asterisk for p-value.
comparison = paste(stage, conditioning, treatment1, "vs", treatment2, sep = ".")
) %>%
select(!t_test_result)
acute_treatment_within_life.stages_conditioning_t_test_results_list[[paste(stage, conditioning, treatment1, treatment2, sep = ".")]] <- t_test_results
}
}
}
}
# Combine results into a single data frame
acute_treatment_within_life.stages_conditioning_t_test_results_df <- bind_rows(acute_treatment_within_life.stages_conditioning_t_test_results_list, .id = "comparison_id")
# View the results
print(acute_treatment_within_life.stages_conditioning_t_test_results_df)
# A tibble: 56 × 6
comparison_id Target estimate_diff p_value asterisk comparison
<chr> <chr> <dbl> <dbl> <chr> <chr>
1 seed.treated.high.ambient ATPsynth… 0.658 0.0602 "" seed.trea…
2 seed.treated.high.ambient DNMT1 -1.55 0.166 "" seed.trea…
3 seed.treated.high.ambient HSP70 -0.133 0.920 "" seed.trea…
4 seed.treated.high.ambient HSP90 0.970 0.0699 "" seed.trea…
5 seed.treated.high.ambient cGAS -0.825 0.0627 "" seed.trea…
6 seed.treated.high.ambient citrate … -1.15 0.0114 "*" seed.trea…
7 seed.treated.high.ambient viperin -0.771 0.255 "" seed.trea…
8 seed.control.ambient.high ATPsynth… -0.304 0.125 "" seed.cont…
9 seed.control.ambient.high DNMT1 -1.54 0.279 "" seed.cont…
10 seed.control.ambient.high HSP70 -1.59 0.0190 "*" seed.cont…
# ℹ 46 more rows
# Create box plots for each comparison
unique_comparisons <- unique(life_stage_t_test_results_df$comparison)
for (comparison in unique_comparisons) {
create_boxplot_delta_Cq(delta_Cq_df, comparison, life_stage_t_test_results_df)
}
# Create box plots for each comparison
unique_comparisons <- unique(conditioning_treatment_t_test_results_df$comparison)
for (comparison in unique_comparisons) {
create_boxplot_conditioning(delta_Cq_df, comparison, conditioning_treatment_t_test_results_df)
}
# Create box plots for each comparison
unique_comparisons <- unique(acute_treatment_t_test_results_df$comparison)
for (comparison in unique_comparisons) {
create_boxplot_acute(delta_Cq_df, comparison, acute_treatment_t_test_results_df)
}
# Loop through each comparison in the t-test results and create box plots
for (comparison in unique(acute_treatment_within_life.stages_conditioning_t_test_results_df$comparison)) {
create_boxplot_acute_conditioning(delta_Cq_df, comparison, acute_treatment_within_life.stages_conditioning_t_test_results_df)
}
# Calculate delta_delta_Cq
delta_delta_conditioning_fold_change <- delta_Cq_df %>%
group_by(life.stage, Target) %>%
summarize(
treated_delta_Cq = mean(delta_Cq[conditioning.treatment == "treated"], na.rm = TRUE),
control_delta_Cq = mean(delta_Cq[conditioning.treatment == "control"], na.rm = TRUE)
) %>%
mutate(delta_delta_Cq = treated_delta_Cq - control_delta_Cq) %>%
select(life.stage, Target, delta_delta_Cq)
str(delta_delta_conditioning_fold_change)
gropd_df [28 × 3] (S3: grouped_df/tbl_df/tbl/data.frame)
$ life.stage : chr [1:28] "adult" "adult" "adult" "adult" ...
$ Target : chr [1:28] "ATPsynthase" "DNMT1" "HSP70" "HSP90" ...
$ delta_delta_Cq: num [1:28] -0.0779 -0.3116 0.941 0.7639 0.1955 ...
- attr(*, "groups")= tibble [4 × 2] (S3: tbl_df/tbl/data.frame)
..$ life.stage: chr [1:4] "adult" "juvenile" "seed" "spat"
..$ .rows : list<int> [1:4]
.. ..$ : int [1:7] 1 2 3 4 5 6 7
.. ..$ : int [1:7] 8 9 10 11 12 13 14
.. ..$ : int [1:7] 15 16 17 18 19 20 21
.. ..$ : int [1:7] 22 23 24 25 26 27 28
.. ..@ ptype: int(0)
..- attr(*, ".drop")= logi TRUE
# Calculate delta_delta_Cq for acute treatment
delta_delta_Cq_acute_df <- delta_Cq_df %>%
group_by(life.stage, Target, acute.treatment) %>%
summarize(
treated_delta_Cq = mean(delta_Cq[conditioning.treatment == "treated"], na.rm = TRUE),
control_delta_Cq = mean(delta_Cq[conditioning.treatment == "control"], na.rm = TRUE)
) %>%
mutate(delta_delta_Cq = treated_delta_Cq - control_delta_Cq) %>%
select(life.stage, Target, acute.treatment, delta_delta_Cq)
str(delta_delta_Cq_acute_df)
gropd_df [56 × 4] (S3: grouped_df/tbl_df/tbl/data.frame)
$ life.stage : chr [1:56] "adult" "adult" "adult" "adult" ...
$ Target : chr [1:56] "ATPsynthase" "ATPsynthase" "DNMT1" "DNMT1" ...
$ acute.treatment: chr [1:56] "ambient" "high" "ambient" "high" ...
$ delta_delta_Cq : num [1:56] -0.112 -0.0438 -0.2467 -0.3765 0.9455 ...
- attr(*, "groups")= tibble [28 × 3] (S3: tbl_df/tbl/data.frame)
..$ life.stage: chr [1:28] "adult" "adult" "adult" "adult" ...
..$ Target : chr [1:28] "ATPsynthase" "DNMT1" "HSP70" "HSP90" ...
..$ .rows : list<int> [1:28]
.. ..$ : int [1:2] 1 2
.. ..$ : int [1:2] 3 4
.. ..$ : int [1:2] 5 6
.. ..$ : int [1:2] 7 8
.. ..$ : int [1:2] 9 10
.. ..$ : int [1:2] 11 12
.. ..$ : int [1:2] 13 14
.. ..$ : int [1:2] 15 16
.. ..$ : int [1:2] 17 18
.. ..$ : int [1:2] 19 20
.. ..$ : int [1:2] 21 22
.. ..$ : int [1:2] 23 24
.. ..$ : int [1:2] 25 26
.. ..$ : int [1:2] 27 28
.. ..$ : int [1:2] 29 30
.. ..$ : int [1:2] 31 32
.. ..$ : int [1:2] 33 34
.. ..$ : int [1:2] 35 36
.. ..$ : int [1:2] 37 38
.. ..$ : int [1:2] 39 40
.. ..$ : int [1:2] 41 42
.. ..$ : int [1:2] 43 44
.. ..$ : int [1:2] 45 46
.. ..$ : int [1:2] 47 48
.. ..$ : int [1:2] 49 50
.. ..$ : int [1:2] 51 52
.. ..$ : int [1:2] 53 54
.. ..$ : int [1:2] 55 56
.. ..@ ptype: int(0)
..- attr(*, ".drop")= logi TRUE
# Calculate delta_delta_Cq for life stage comparisons
delta_delta_Cq_life_stage_df <- delta_Cq_df %>%
group_by(Target, life.stage) %>%
summarize(mean_delta_Cq = mean(delta_Cq, na.rm = TRUE)) %>%
ungroup() %>%
pivot_wider(names_from = life.stage, values_from = mean_delta_Cq) %>%
mutate(
delta_delta_Cq_adult_vs_seed = adult - seed,
delta_delta_Cq_spat_vs_seed = spat - seed,
delta_delta_Cq_adult_vs_spat = adult - spat
) %>%
pivot_longer(cols = starts_with("delta_delta_Cq_"), names_to = "comparison", values_to = "delta_delta_Cq") %>%
filter(!is.na(delta_delta_Cq))
# Display the structure of the resulting data frame
str(delta_delta_Cq_life_stage_df)
tibble [21 × 7] (S3: tbl_df/tbl/data.frame)
$ Target : chr [1:21] "ATPsynthase" "ATPsynthase" "ATPsynthase" "DNMT1" ...
$ adult : num [1:21] 0.48 0.48 0.48 6.17 6.17 ...
$ juvenile : num [1:21] 0.378 0.378 0.378 5.808 5.808 ...
$ seed : num [1:21] 0.17 0.17 0.17 6.24 6.24 ...
$ spat : num [1:21] 0.337 0.337 0.337 6.443 6.443 ...
$ comparison : chr [1:21] "delta_delta_Cq_adult_vs_seed" "delta_delta_Cq_spat_vs_seed" "delta_delta_Cq_adult_vs_spat" "delta_delta_Cq_adult_vs_seed" ...
$ delta_delta_Cq: num [1:21] 0.3098 0.1664 0.1433 -0.0736 0.2019 ...
# Calculate delta_delta_Cq for acute treatment comparisons within each life stage and conditioning treatment
delta_delta_Cq_acute_within_life_stage_conditioning_df <- delta_Cq_df %>%
group_by(life.stage, conditioning.treatment, Target, acute.treatment) %>%
summarize(mean_delta_Cq = mean(delta_Cq, na.rm = TRUE)) %>%
ungroup() %>%
pivot_wider(names_from = acute.treatment, values_from = mean_delta_Cq) %>%
mutate(delta_delta_Cq_high_vs_ambient = high - ambient) %>%
pivot_longer(cols = starts_with("delta_delta_Cq_"), names_to = "comparison", values_to = "delta_delta_Cq") %>%
filter(!is.na(delta_delta_Cq))
# Display the structure of the resulting data frame
str(delta_delta_Cq_acute_within_life_stage_conditioning_df)
tibble [56 × 7] (S3: tbl_df/tbl/data.frame)
$ life.stage : chr [1:56] "adult" "adult" "adult" "adult" ...
$ conditioning.treatment: chr [1:56] "control" "control" "control" "control" ...
$ Target : chr [1:56] "ATPsynthase" "DNMT1" "HSP70" "HSP90" ...
$ ambient : num [1:56] 0.566 6.448 3.944 1.259 5.207 ...
$ high : num [1:56] 0.472 6.199 3.673 0.29 4.609 ...
$ comparison : chr [1:56] "delta_delta_Cq_high_vs_ambient" "delta_delta_Cq_high_vs_ambient" "delta_delta_Cq_high_vs_ambient" "delta_delta_Cq_high_vs_ambient" ...
$ delta_delta_Cq : num [1:56] -0.0946 -0.2492 -0.2715 -0.969 -0.5983 ...
# Calculate fold change and output to a new data frame
fold_change_life_stage_df <- delta_delta_Cq_life_stage_df %>%
mutate(fold_change = 2^(-delta_delta_Cq))
# Display the structure of the resulting data frame
str(fold_change_life_stage_df)
tibble [21 × 8] (S3: tbl_df/tbl/data.frame)
$ Target : chr [1:21] "ATPsynthase" "ATPsynthase" "ATPsynthase" "DNMT1" ...
$ adult : num [1:21] 0.48 0.48 0.48 6.17 6.17 ...
$ juvenile : num [1:21] 0.378 0.378 0.378 5.808 5.808 ...
$ seed : num [1:21] 0.17 0.17 0.17 6.24 6.24 ...
$ spat : num [1:21] 0.337 0.337 0.337 6.443 6.443 ...
$ comparison : chr [1:21] "delta_delta_Cq_adult_vs_seed" "delta_delta_Cq_spat_vs_seed" "delta_delta_Cq_adult_vs_spat" "delta_delta_Cq_adult_vs_seed" ...
$ delta_delta_Cq: num [1:21] 0.3098 0.1664 0.1433 -0.0736 0.2019 ...
$ fold_change : num [1:21] 0.807 0.891 0.905 1.052 0.869 ...
delta_delta_conditioning_fold_change <- delta_delta_conditioning_fold_change %>%
mutate(fold_change = 2^(-delta_delta_Cq)) %>%
distinct(Target, fold_change)
str(delta_delta_conditioning_fold_change)
gropd_df [28 × 3] (S3: grouped_df/tbl_df/tbl/data.frame)
$ life.stage : chr [1:28] "adult" "adult" "adult" "adult" ...
$ Target : chr [1:28] "ATPsynthase" "DNMT1" "HSP70" "HSP90" ...
$ fold_change: num [1:28] 1.055 1.241 0.521 0.589 0.873 ...
- attr(*, "groups")= tibble [4 × 2] (S3: tbl_df/tbl/data.frame)
..$ life.stage: chr [1:4] "adult" "juvenile" "seed" "spat"
..$ .rows : list<int> [1:4]
.. ..$ : int [1:7] 1 2 3 4 5 6 7
.. ..$ : int [1:7] 8 9 10 11 12 13 14
.. ..$ : int [1:7] 15 16 17 18 19 20 21
.. ..$ : int [1:7] 22 23 24 25 26 27 28
.. ..@ ptype: int(0)
..- attr(*, ".drop")= logi TRUE
# Calculate fold change for acute treatment
delta_delta_acute_fold_change <- delta_delta_Cq_acute_df %>%
mutate(fold_change = 2^(-delta_delta_Cq)) %>%
distinct(life.stage, Target, acute.treatment, fold_change)
# Display the structure of the resulting data frame
str(delta_delta_acute_fold_change)
gropd_df [56 × 4] (S3: grouped_df/tbl_df/tbl/data.frame)
$ life.stage : chr [1:56] "adult" "adult" "adult" "adult" ...
$ Target : chr [1:56] "ATPsynthase" "ATPsynthase" "DNMT1" "DNMT1" ...
$ acute.treatment: chr [1:56] "ambient" "high" "ambient" "high" ...
$ fold_change : num [1:56] 1.081 1.031 1.186 1.298 0.519 ...
- attr(*, "groups")= tibble [28 × 3] (S3: tbl_df/tbl/data.frame)
..$ life.stage: chr [1:28] "adult" "adult" "adult" "adult" ...
..$ Target : chr [1:28] "ATPsynthase" "DNMT1" "HSP70" "HSP90" ...
..$ .rows : list<int> [1:28]
.. ..$ : int [1:2] 1 2
.. ..$ : int [1:2] 3 4
.. ..$ : int [1:2] 5 6
.. ..$ : int [1:2] 7 8
.. ..$ : int [1:2] 9 10
.. ..$ : int [1:2] 11 12
.. ..$ : int [1:2] 13 14
.. ..$ : int [1:2] 15 16
.. ..$ : int [1:2] 17 18
.. ..$ : int [1:2] 19 20
.. ..$ : int [1:2] 21 22
.. ..$ : int [1:2] 23 24
.. ..$ : int [1:2] 25 26
.. ..$ : int [1:2] 27 28
.. ..$ : int [1:2] 29 30
.. ..$ : int [1:2] 31 32
.. ..$ : int [1:2] 33 34
.. ..$ : int [1:2] 35 36
.. ..$ : int [1:2] 37 38
.. ..$ : int [1:2] 39 40
.. ..$ : int [1:2] 41 42
.. ..$ : int [1:2] 43 44
.. ..$ : int [1:2] 45 46
.. ..$ : int [1:2] 47 48
.. ..$ : int [1:2] 49 50
.. ..$ : int [1:2] 51 52
.. ..$ : int [1:2] 53 54
.. ..$ : int [1:2] 55 56
.. ..@ ptype: int(0)
..- attr(*, ".drop")= logi TRUE
# Calculate fold change for acute treatment comparisons within each life stage and conditioning treatment
fold_change_acute_within_life_stage_conditioning_df <- delta_delta_Cq_acute_within_life_stage_conditioning_df %>%
mutate(fold_change = 2^(-delta_delta_Cq))
# Display the structure of the resulting data frame
str(fold_change_acute_within_life_stage_conditioning_df)
tibble [56 × 8] (S3: tbl_df/tbl/data.frame)
$ life.stage : chr [1:56] "adult" "adult" "adult" "adult" ...
$ conditioning.treatment: chr [1:56] "control" "control" "control" "control" ...
$ Target : chr [1:56] "ATPsynthase" "DNMT1" "HSP70" "HSP90" ...
$ ambient : num [1:56] 0.566 6.448 3.944 1.259 5.207 ...
$ high : num [1:56] 0.472 6.199 3.673 0.29 4.609 ...
$ comparison : chr [1:56] "delta_delta_Cq_high_vs_ambient" "delta_delta_Cq_high_vs_ambient" "delta_delta_Cq_high_vs_ambient" "delta_delta_Cq_high_vs_ambient" ...
$ delta_delta_Cq : num [1:56] -0.0946 -0.2492 -0.2715 -0.969 -0.5983 ...
$ fold_change : num [1:56] 1.07 1.19 1.21 1.96 1.51 ...
library(ggplot2)
# Generate bar plots for each group of comparison within each life stage and conditioning treatment
plot_list <- fold_change_acute_within_life_stage_conditioning_df %>%
split(list(.$life.stage, .$conditioning.treatment, .$comparison)) %>%
lapply(function(df) {
life_stage <- unique(df$life.stage)
conditioning_treatment <- unique(df$conditioning.treatment)
comparison_title <- gsub("delta_delta_Cq_", "", unique(df$comparison))
comparison_title <- gsub("_vs_", " vs. ", comparison_title)
ggplot(df, aes(x = Target, y = fold_change)) +
geom_bar(stat = "identity") +
labs(title = paste("Gene Expression -", life_stage, "-", conditioning_treatment, "-", comparison_title),
x = "Target", y = "Fold Change") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
})
# Display the plots
for (plot in plot_list) {
print(plot)
}
# Generate bar plots for each group of comparison
plot_list <- fold_change_life_stage_df %>%
split(.$comparison) %>%
lapply(function(df) {
comparison_title <- gsub("delta_delta_Cq_", "", unique(df$comparison))
comparison_title <- gsub("_vs_", " vs. ", comparison_title)
ggplot(df, aes(x = Target, y = fold_change)) +
geom_bar(stat = "identity") +
labs(title = paste("Gene Expression -", comparison_title), x = "Target", y = "Fold Change") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
})
# Display the plots
for (plot in plot_list) {
print(plot)
}
hist(delta_Cq_df$delta_Cq)
Run an anova model to test for effects of lifestage, conditioning, and acute treatment on delta Cq values for each target.
ATP synthase
library(car)
model<-delta_Cq_df%>%
filter(Target=="ATPsynthase")%>%
aov(delta_Cq ~ life.stage * conditioning.treatment * acute.treatment, data=.)
summary(model)
qqPlot(model$residuals)
No significant effects.
DNMT1
model<-delta_Cq_df%>%
filter(Target=="DNMT1")%>%
aov(delta_Cq ~ life.stage * conditioning.treatment * acute.treatment, data=.)
summary(model)
qqPlot(model$residuals)
No significant effects.
HSP70
model<-delta_Cq_df%>%
filter(Target=="HSP70")%>%
aov(delta_Cq ~ life.stage * conditioning.treatment * acute.treatment, data=.)
summary(model)
qqPlot(model$residuals)
Significant effect of conditioning x acute treatment.
HSP90
model<-delta_Cq_df%>%
filter(Target=="HSP90")%>%
aov(delta_Cq ~ life.stage * conditioning.treatment * acute.treatment, data=.)
summary(model)
qqPlot(model$residuals)
Significant effect of lifestage x acute treatment, lifestage x conditioning treatment, acute treatment, and lifestage.
cGAS
model<-delta_Cq_df%>%
filter(Target=="cGAS")%>%
aov(delta_Cq ~ life.stage * conditioning.treatment * acute.treatment, data=.)
summary(model)
qqPlot(model$residuals)
Significant effect of lifestage x acute treatment and lifestage.
Citrate synthase
model<-delta_Cq_df%>%
filter(Target=="citrate synthase")%>%
aov(delta_Cq ~ life.stage * conditioning.treatment * acute.treatment, data=.)
summary(model)
qqPlot(model$residuals)
Significant effect of acute treatment and lifestage x acute treatment.
Display plot of acute x conditioning treatment faceted for each lifestage for each target.
plot1<-delta_Cq_df%>%
filter(Target=="ATPsynthase")%>%
group_by(Target, life.stage, acute.treatment, conditioning.treatment)%>%
summarise(mean=mean(delta_Cq, na.rm=TRUE), se=sd(delta_Cq, na.rm=TRUE)/sqrt(length(delta_Cq)))%>%
ggplot(aes(x=acute.treatment, y=mean, colour=conditioning.treatment))+
facet_grid(~life.stage)+
scale_colour_manual(values=c("darkgray", "orange"))+
geom_point()+
geom_line(aes(group=conditioning.treatment))+
geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=0.1)+
ggtitle("ATP synthase (no sign. effects)")+
ylab("Delta Cq")+
ylim(-0.5,1.5)+
geom_hline(yintercept=0, linetype="dashed")+
theme_classic();plot1
plot2<-delta_Cq_df%>%
filter(Target=="DNMT1")%>%
group_by(Target, life.stage, acute.treatment, conditioning.treatment)%>%
summarise(mean=mean(delta_Cq, na.rm=TRUE), se=sd(delta_Cq, na.rm=TRUE)/sqrt(length(delta_Cq)))%>%
ggplot(aes(x=acute.treatment, y=mean, colour=conditioning.treatment))+
facet_grid(~life.stage)+
geom_point()+
scale_colour_manual(values=c("darkgray", "orange"))+
geom_line(aes(group=conditioning.treatment))+
geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=0.1)+
ggtitle("DNMT (no sign. effects)")+
ylab("Delta Cq")+
geom_hline(yintercept=0, linetype="dashed")+
ylim(0,9)+
theme_classic();plot2
plot3<-delta_Cq_df%>%
filter(Target=="HSP70")%>%
group_by(Target, life.stage, acute.treatment, conditioning.treatment)%>%
summarise(mean=mean(delta_Cq, na.rm=TRUE), se=sd(delta_Cq, na.rm=TRUE)/sqrt(length(delta_Cq)))%>%
ggplot(aes(x=acute.treatment, y=mean, colour=conditioning.treatment))+
facet_grid(~life.stage)+
scale_colour_manual(values=c("darkgray", "orange"))+
geom_point()+
geom_line(aes(group=conditioning.treatment))+
geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=0.1)+
ggtitle("HSP70 (sign. conditioning x acute treatment)")+
ylab("Delta Cq")+
ylim(0,8)+
geom_hline(yintercept=0, linetype="dashed")+
theme_classic();plot3
plot3a<-delta_Cq_df%>%
filter(Target=="HSP70")%>%
group_by(Target, acute.treatment, conditioning.treatment)%>%
summarise(mean=mean(delta_Cq, na.rm=TRUE), se=sd(delta_Cq, na.rm=TRUE)/sqrt(length(delta_Cq)))%>%
ggplot(aes(x=acute.treatment, y=mean, colour=conditioning.treatment))+
geom_point()+
scale_colour_manual(values=c("darkgray", "orange"))+
geom_line(aes(group=conditioning.treatment))+
geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=0.1)+
ggtitle("HSP70 (sign. conditioning x acute treatment)")+
ylab("Delta Cq")+
geom_hline(yintercept=0, linetype="dashed")+
ylim(0,8)+
theme_classic();plot3a
plot4<-delta_Cq_df%>%
filter(Target=="HSP90")%>%
group_by(Target, life.stage, acute.treatment, conditioning.treatment)%>%
summarise(mean=mean(delta_Cq, na.rm=TRUE), se=sd(delta_Cq, na.rm=TRUE)/sqrt(length(delta_Cq)))%>%
ggplot(aes(x=acute.treatment, y=mean, colour=conditioning.treatment))+
facet_grid(~life.stage)+
scale_colour_manual(values=c("darkgray", "orange"))+
geom_point()+
geom_line(aes(group=conditioning.treatment))+
geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=0.1)+
ggtitle("HSP90 (sign. lifestage x acute & lifestage x conditioning treatment)")+
ylab("Delta Cq")+
ylim(-2,2.5)+
geom_hline(yintercept=0, linetype="dashed")+
theme_classic();plot4
plot5<-delta_Cq_df%>%
filter(Target=="cGAS")%>%
group_by(Target, life.stage, acute.treatment, conditioning.treatment)%>%
summarise(mean=mean(delta_Cq, na.rm=TRUE), se=sd(delta_Cq, na.rm=TRUE)/sqrt(length(delta_Cq)))%>%
ggplot(aes(x=acute.treatment, y=mean, colour=conditioning.treatment))+
facet_grid(~life.stage)+
scale_colour_manual(values=c("darkgray", "orange"))+
geom_point()+
geom_line(aes(group=conditioning.treatment))+
geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=0.1)+
ggtitle("cGAS (sign. lifestage x acute treatment)")+
ylab("Delta Cq")+
ylim(2,8)+
geom_hline(yintercept=0, linetype="dashed")+
theme_classic();plot5
plot6<-delta_Cq_df%>%
filter(Target=="citrate synthase")%>%
group_by(Target, life.stage, acute.treatment, conditioning.treatment)%>%
summarise(mean=mean(delta_Cq, na.rm=TRUE), se=sd(delta_Cq, na.rm=TRUE)/sqrt(length(delta_Cq)))%>%
ggplot(aes(x=acute.treatment, y=mean, colour=conditioning.treatment))+
facet_grid(~life.stage)+
scale_colour_manual(values=c("darkgray", "orange"))+
geom_point()+
geom_line(aes(group=conditioning.treatment))+
geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=0.1)+
ggtitle("cGAS (sign. lifestage x acute treatment)")+
ylab("Delta Cq")+
ylim(-2,3)+
geom_hline(yintercept=0, linetype="dashed")+
theme_classic();plot6
The HSP genes show effects of conditioning. The other genes show effects of lifestage and/or acute stress.
Lifestage specific stress responses.